#include "density_matrix.hpp"


template<typename T>
QuantumWalk::density_matrix<T>::density_matrix(const state<T>& psi):
  gen_op(psi.no_vertices(),psi.no_internal_deg()){
  //Actual initialization
  int N=this->dimV_ * this->dimD_;
  data_.resize(N,N);
  //rho(k,l)=rho(ij,vw), where k=i*dimV_+v  and l=j*dimV_+w
  for(int i=0;i<this->dimD_;i++){//Iterate internal degrees
    for(int j=0;j<this->dimD_;j++){//Iterate internal degrees
      for(int v=0;v<this->dimV_;v++){//Iterate sites
	for(int w=0;w<this->dimV_;w++){//iterate sites
	  data_(i*this->dimV_+v,j*this->dimV_+w)=
	    std::conj(psi(i,v))*psi(j,w);
	}
      }
    }
  }
};

template<typename T> 
QuantumWalk::density_matrix<T>::density_matrix(matrix<T> m):
  gen_op(m.size1(),m.size2()){
  int N=this->dimV_ * this->dimD_;
  data_.resize(N,N);
  for(int i=0;i<this->dimD_;i++){//Iterate internal degrees
    for(int j=0;j<this->dimD_;j++){//Iterate internal degrees
      for(int v=0;v<this->dimV_;v++){//Iterate sites
	for(int w=0;w<this->dimV_;w++){//iterate sites
	  int irow,icol;
	  irow=this->two_to_one(i,v);
	  icol=this->two_to_one(j,w);
	  data_(icol,irow)=
	    std::conj(m(i,v))*m(j,w);
	}
      }
    }
  }
};

template<typename T>
matrix<T>& QuantumWalk::density_matrix<T>::asmatrix(){
  return data_;
};

template class QuantumWalk::density_matrix<std::complex<double> >;
